Kinetically driven ordered phase formation in binary colloidal crystals 
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The aggregation of binary colloids of same size and balanced charges is studied by Brownian dynamics 
simulations for dilute suspensions. It is shown that, under appropriate conditions, the formation of colloidal 
crystals is dominated by kinetic effects leading to the growth of well-ordered crystallites of the sodium-chloride 
(NaCl) bulk phase. These crystallites form with very high probability even when the cesium-chloride (CsCl) 
phase is more stable thermodynamically. Global optimization searches show that this result is not related to the 
most favorable structures of small clusters, that are either amorphous or of CsCl structure. The formation of 
the NaCl phase is related to the specific kinetics of the crystallization process, which takes place by a two-step 
mechanism. In this mechanism, dense fluid aggregates form at first and then crystallization follows. It is shown 
that the type of short-range order in these dense fluid aggregates determines which phase is finally formed in 
the crystallites. The role of hydrodynamic effects in the aggregation process is analyzed by Stochastic Rotation 
Dynamics - Molecular Dynamics simulations, finding that these effects do not play a major role in the formation 
of the crystallites. 



INTRODUCTION 

Colloidal suspensions have been widely studied for a vari- 
ety of applications. From the point of view of basic science, 
a growing interest stems from the fact that colloidal suspen- 
sions are many-particle systems that are easier to observe than 
molecular systems. Aggregation and nucleation phenomena 
of colloids are therefore simpler to analyze, also by means 
of real-space techniques. Moreover, the interactions between 
colloids can be tuned in a controlled way so that a variety of 
structures with different properties can be experimentally pro- 
duced [lfj. Recently, much attention has been devoted to bi- 
nary colloids in which two types of colloidal particles (types 
A and B in the following) can acquire opposite charges in 
suspension, interacting mainly through screened electrostatic 
forces. Binary colloids can thus aggregate J^j-Ql, in a pro- 
cess which is usually called heteroaggregation. The compe- 
tition between attractive and repulsive forces gives rise to a 
rich phenomenology, which is the subject of active investiga- 
tion. In fact, binary colloids have been shown to be able to 
form either compact colloidal crystals with a good degree of 
long-range order HI El H> or colloidalgels characterized by a 
network of fractal-like aggregates []7j-|9|],the latter being more 
tenuous than in single-component gels 19(]- 

From the theoretical/computational point of view, notable 
attention has been devoted both to the determination of the 
equilibrium properties of binary suspensions J2I Q [lQ 11 1 . 
with the calculation of phase diagrams of representative model 
systems, and to their aggregation kinetics J2, EH-G3] ■ How- 
ever, several aspects of heteroaggregation are still to be un- 
derstood, especially for what concerns the basic mechanisms 
that determine the structures obtained in the actual aggrega- 
tion processes. 

An important issue in this colloidal aggregation is to deter- 



mine to what extent the structures that are actually formed in 
the growth of colloidal crystals reflect their thermodynamic 
equilibrium structures. This has been investigated in the 
case of the heteroaggregation of colloids interacting through 
screened electrostatic forces by Sanz et al. IU2II . in the case of 
high volume fraction $ ~ 0.47. In their simulations, Sanz et 
al. used the forward-flux sampling method Ill5ll to determine 
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nucleation rates. They found that highest nucleation rate was 
associated to the formation of a substitutionally disordered fee 
crystal, which is neither the lowest free-energy phase nor the 
phase characterized by the lowest free-energy barrier for nu- 
cleation from the melt. In fact, the lowest free-energy struc- 
ture is of cesium-chloride (CsCl) type. Sanz et al. attributed 
this finding to the fact that sub-critical small clusters have dis- 
ordered fee structures, so that growth continues keeping this 
kind of structure because the transition to the CsCl structure 
is kinetically inhibited. 

In this article we perform Brownian dynamics simulations 
of crystal nucleation and growth in the same model system as 
Sanz et. al. consisting of binary spherical colloids 
of same radius a and absolute value of the charge, interact- 
ing through a screened electrostatic potential of Yukawa-like 
form. This potential gives a valid description of the interac- 
tions between colloids in index-matched solvents IH floll . For 
a given colloid volume fraction <£>, the equilibrium behavior 
of colloids interacting by this potential can be described as a 
function of two parameters: U* = Uq/ {ksT), where Uq is the 
well depth, and kcl, where k _1 is the Debye screening length. 
At variance with Refs. HI 0, we consider a quite dilute 
system, with a volume fraction $ = 0.05, which is below the 
percolation threshold. 

Our simulations show that aggregates corresponding to 
metastable colloidal crystals form with very high probability 
under appropriate conditions. While in the case of high vol- 
ume fractions treated in Ref. IU2I1 the metastable phase was 
disordered, in our dilute system we find that a very well or- 
dered metastable phase form. Specifically, we show that crys- 
tallites with the sodium-chloride (NaCl) structure grow for 



2 



values of (U* , na) at which the cesium-chloride (CsCl) crys- 
tal structure is the most stable from the thermodynamic point 
of view. The formation of the NaCl-type crystallites is ratio- 
nalized in the framework of the two-step mechanism of crystal 
nucleation | ID- 1.8], in which a dense colloidal fluid acts as a 
precursor for crystal nucleation. We show specifically that the 
type of short-range order in this dense fluid is crucial for de- 
termining which crystal phase is going to be formed. We show 
also that the NaCl phase formation is not related to energet- 
ically favorable structures of small clusters, the latter being 
either amorphous or with CsCl structure. Colloidal crystals 
with the NaCl structure were obtained with colloids of differ- 
ent sizes Il9ll . Our results demonstrate these crystals can grow 
in systems with balanced size and charge. We analyze also the 
importance of hydrodynamic effects in the nucleation of the 
NaCl phase by performing Stochastic Rotation Dynamics - 
Molecular Dynamics (SRD-MD) simulations which 
show that these effects do not play a major role. 



MODEL AND METHODS 

In our simulations the interaction potential between the col- 
loids is of the form 



2a 

Ui 3 {r) = U sign(qiqj)—e 



-«(r-2o) 



(1) 



where Uq is a constant, which gives the value of the potential 
at contact, and the sign of the charge qi depends on the col- 
loid i and being either of type A or B. In the following we 
use U* = Uo/(kBT) = 9. This expression of the potential is 
valid for r > 2a. For r < 2a, the hard-wall repulsion is ap- 
proximated by a power law increasing as We have 
checked that, in infinite crystals, the Madelung energy favors 
the CsCl phase over the NaCl phase for all values of na, as 
pointed out in Ref. U]. The Madelung energy difference be- 
tween the two phases decreases with decreasing na. In Ref. 
|l|] it was also shown that, at finite temperature, entropic con- 
tributions increase the stability of the NaCl bulk phase, which 
should however become more stable than the CsCl phase only 
below na ~ 2.5 for U* of the magnitude considered here. 

In the following, we use three different simulation tech- 
niques, Brownian Dynamics (BD) simulations, Basin Hop- 
ping (BH) global optimization searches, and SRD-MD sim- 
ulations. 

Our BD simulations are performed by numerically solv- 
ing the Langevin equation by means of the algorithm de- 
scribed in Ref. (Ht]. We choose the dimensionless time unit as 
t* = 4a 2 /Dq, where Dq is the diffusion coefficient of isolated 
colloids in the suspension, t* is the time taken by an isolated 
colloid to diffuse across a distance equivalent to its diame- 
ter. Do is given by Einstein's relation, in which the solvent 
friction and the colloid mass are chosen as in Ref. [@]. The 
Langevin equation is solved by using time steps St in the range 
10~ 5 — 10~ 6 t*. The number of colloids N in the simulations 
is 200, half being positively and half negatively charged. The 



colloids are initially placed in a cubic box of volume V with 
random positions. Periodic boundary conditions are applied. 
Simulations are stopped after a time i/j„ ~ 2 10 3 i*, which 
is sufficient to let all colloids form a single aggregate of crys- 
talline structure. For comparison, in the system of Ref. fl8[, 
this tf i n would amount to ~ 10 3 s. 

The lowest-energy structures of colloidal clusters have been 
searched for by a global optimization procedure employing 
the BH method 12311 . which has proven to be very efficient 
in the global optimization of atomic clusters. This method 
consists of a Metropolis Monte Carlo procedure in which the 
energy is locally minimized at each step. BH has been ap- 
plied recently also to colloidal clusters in Ref. Il24ll, w here 
the computational procedure is described (see also II25II for 
a more complete description of the method). We note here 
that we have used three different kinds of elementary moves 
to generate new configurations in the Monte Carlo procedure. 
These are the Brownian move, the shake move |25l|, and a 
combined move, made of an uniform compression, in which 
the whole cluster is compressed by a factor 0.8 - 0.9, followed 
by a mild shake, i.e. of small random displacements of all 
colloids around their positions. After each move, the cluster 
is locally relaxed to generate the new configuration, which is 
then either accepted or refused according to the Metropolis 
rule. In all cases, 5 simulations of 10 5 steps each have been 
performed. 

Hydrodynamic effects have been taken into account by 
means of SRD-MD simulations. SRD-MD is a combination 
of Stochastic Rotation Dynamics and standard Molecular Dy- 
namics that allows taking into account long range hydrody- 
namic effects in the study of colloidal suspensions H201 12111 . 
This technique has been chosen for its relative efficiency and 
low computational cost, to check whether hydrodynamic ef- 
fects could affect the nucleation process in our system. SRD- 
MD is based on the use of fake fluid particles that mimic the 
fluid properties (essentially its momentum). The simulation 
box is divided into smaller cubic cells. In each cell the fluid 
particles exchange momentum thanks to a rotation of their ve- 
locities. For an accurate description of the method we refer to 
the literature |2o|-22l]- Here we report only the choices of the 
most important parameters for our system. The linear size of 
the cells ao is taken as half of the radius of the colloids. For 
the key parameters of the SRD-MD simulations we made the 
same choice proposed by Padding and Louis [21]: the num- 
ber of fluid particles per cell 7 equal to 5, the dimensionless 
mean free path A equal to 0. 1 and the rotation angle a equal 
to 90°. The parameters have been chosen to reproduce the 
same diffusion coefficient of isolated colloids in both BD and 



SRD-MD Il22h . Recently it has been proved that these choices 
are able to correctly reproduce the essential features of the 
hydrodynamic interactions in a colloidal suspension II22I1 . For 
what concerns the energetic model, we have chosen the same 
Yukawa potential used in BD simulations for colloid-colloid 
interactions, while for colloid-fluid interactions we have cho- 
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sen a soft-sphere potential 

V cf (r) = £ cf 
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for r < r c = 2.5 <r c f and otherwise. We have chosen n=l2 
and e c f = 2.5 ksT. The value of a c f has been taken smaller 
than the colloidal radius (a c f = 0.8 a) to avoid spurious de- 
pletion attractions between colloids |21 ] . 



RESULTS AND DISCUSSION 
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TABLE I: Final structure of the aggregates obtained in Brownian 
simulations for different values of na and volume fraction $. Either 
10 or 20 simulations have been performed in each case. 



Let us consider the aggregation of the colloids for na = 3, 
at which the NaCl crystal phase is never found to be the most 
stable [Ql QjJ from the thermodynamic point of view. The 
results concerning final aggregates in BD simulations (see Ta- 
ble HI show that in a large majority of cases (15 in 20 simu- 
lations) structures of NaCl type are formed. The remaining 
simulations give either mixed or CsCl structures, the latter 
being found only in two cases. Representative snapshots of 
NaCl, mixed and CsCl structures are shown in Fig. Q] Also 
the simulations for $ = 0.08 give a large majority of NaCl- 
type aggregates. As Ka increases the CsCl phase becomes 
more and more energetically favorable. Therefore, we expect 
that there should be a critical value (na) c at which the proba- 
bility of forming both types of aggregates is the same. From 
the results in Table H] it turns out that (k«) c is close to 3.3, 
a value which is significantly larger than the Ka at which the 
bulk NaCl and CsCl phases coexist in thermodynamic equi- 
librium [[it]. Increasing na further leads to the formation of 
CsCl-type aggregates with an overwhelming probability. On 
the other hand, if Ka is decreased to 2.55, NaCl-type crystal- 
lites are always formed. 



for the lowest-energy structures of colloidal clusters in the size 
range from A 8 B 8 to A32B32. The main structural motifs sin- 
gled out for the Ka values reported in Tableware shown in Fig. 
12 where a few selected sizes are considered (AgBg, A9B9, 
A10B10, AnBii, A13B13, A20B20, and A32B32). Depending 
on Ka, the dominant motif may change as follows. 






FIG. 1: . From left to right, aggregates of NaCl, mixed and CsCl 
structures. Particles of different charges are shown in green (dark 
grey) and red (light grey). 

What are the driving forces leading to the formation of these 
aggregates of the metastable NaCl bulk phase? 

A possible explanation of this result may originate from the 
fact that the aggregation process starts with the formation of 
very small colloidal clusters that subsequently grow and coa- 
lesce, as shown in the simulation of the aggregation in a sim- 
ilar system . Therefore, the structure of the final aggreg ate 
may preserve memory of the equilibrium structures of small 
aggregates. In order to investigate this point we have searched 



FIG. 2: Lowest-energy structures of the aggregates as determined 
by global optimiziation searches, in the cases of AgBa (al and a2), 
A9B9 (bl and b2), Ai Bi (cl and c2), A11B11 (dl and d2), A13B13 
(el and e2), A20B20 (f) and A32B32 (gl and g2). The actual lowest- 
energy structure depends on Ka as explained in the text. 

For small aggregates, from AgBs to A10B10, the lowest- 
energy structures either belong to the CsCl motif, as struc- 
tures al, bl, and cl in Fig. [2] or are disordered (i.e. they are 
neither of CsCl- nor of NaCl-type, see a2, b2, c2 in Fig. [2). 
CsCl-type structures prevail at Ka > 5, while at smaller Ka 
disordered structures are more favourable. From AnBn on, 
CsCl-type structures (dl, el, f, gl in Fig. |2]i always prevail 
with a single exception which is found at Ka = 2.55 only. For 
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this kcl, the best structure of A13B13 is of NaCl type, being 
a 3x3 cube with one vertex missing (structure e2). However, 
with increasing cluster size NaCl structures become more and 
more unfavorable. In fact, at the next magic size for NaCl 
structures, i.e. at size 64, where a perfect NaCl cube A32B32 
can be formed (structure gl), the CsCl structure is lower in 
energy for all ko, down to 2.55. 

These results show that structures of CsCl type prevail al- 
ready for small aggregates, indicating that there are no evi- 
dent effects that favor the NaCl structures at small sizes with 
respect to what is found in bulk crystals. Therefore the expla- 
nation of the growth of NaCl-type aggregates originates from 
a different effect, which is better understood by analyzing in 
details the aggregation process in the simulations. 

To this end, we define two order parameters that are able to 
discriminate between the different crystalline phases. These 
parameters (referred to as Pi and P2 in the following) are de- 
fined as follows. Pi is the fraction of colloids having more 
than 6 first neighbors. This indicates the possible presence 
of CsCl-type short range order, because in NaCl-type aggre- 
gates the number of first neighbors of a given colloid cannot 
be larger than 6. On the other hand, P2 is related to the con- 
figuration of second neighbors, being defined as 
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where is the distance between colloids i and j (i and j 
being of the same type), r2N = 2y/2a and T2c = 4a/-\/3 
are the average equilibrium distance of second neighbors in 
the NaCl and CsCl structures, and a takes into account the 
effect of vibrations (in our simulations we choose a = 0.15a, 
a value that is more than three times smaller than r2N — T2c)- 
Positive and negative values of P2 respectively indicate the 
prevalence of either NaCl- or CsCl-type short-range order in 
the aggregate. In perfect NaCl and CsCl infinite crystals P2 is 
close to 1 and -1, respectively. 

Let us consider the typical behavior of a simulation for 
na = 3 (see Fig. [3] crosses). After a rather short time, 
essentially all colloids attach to some aggregate and, as in 
most cases, a single aggregate is quickly formed (typically for 
t < 100t*). Pi stays very close to 0, showing that the prob- 
ability of configurations in which a colloid has more than six 
first neighbours is very low. On the other hand, P2 increases 
from the initial value of to a positive value P2 ~ 0.05, indi- 
cating that some NaCl-type short-range order is present. How- 
ever, if we look at the aggregate itself, we find configurations 
as those reported in the snapshot al of Fig. |4] The aggregate 
does not show any long-range order, but a liquid-like struc- 
ture, in which some colloids surrounded by a NaCl-type short- 
range order can be singled out. After a lengthy evolution in 
this disordered phase, there is a rather sudden increase of P2, 
which jumps to a value of 0.4, corresponding to the transfor- 
mation of the disordered aggregate into an ordered NaCl-type 
crystallite, as the one shown in snapshot a2 of Fig. 4. 

This behavior can be compared to the case of na = 5. 
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FIG. 3: Behavior of order parameters Pi and P2 as functions of time 
for representative simulations at different ko: 2.55 (circles), 3.00 
(crosses), 5.00 (triangles) and 9.00 (squares). Time is measured in 
units of t* . The curves are related to single simulations, but both Pi 
and P2 have been averaged each 10 3 steps in order to give smoother 
behaviors. 



Again, a single aggregate is rather quickly formed. The ag- 
gregate does not show long-range order, but a somewhat pre- 
vailing CsCl-type short-range order (see bl in Fig. |4} cor- 
responding to a negative P2, and to a small but significantly 
positive Pi . After some time, the sudden transition to the or- 
dered crystallite takes place, but this time the structure is of 
CsCl type (b2 in Fig. g}. 

This analysis shows that the final structure of the crystal- 
lite is therefore well described by the two-step mechanism of 
crystal nucleation |3 12 I- Starting from diluted configura- 
tions ($ < 0.08) dense aggregates form with a metastable 
liquid-like structure which survives for some time. However, 
a transition to an ordered structure finally takes place, and the 
type of ordered structure is correlated, with very high prob a- 
bility, with the features of short-range order in the metastable 
liquid, as described for example by the parameter P2- Specif- 
ically we find that the final structure is almost always of NaCl 
type for disordered aggregates in which P2 > 0.05, while for 
P-2 < 0.02 the CsCl structure is found in most cases. This 



means that the final ordered structure can be predicted with 
high probability by inspecting the relative positions of the col- 
loids in the disordered state. We note that this is not always 
the case in colloid crystallization, as demonstrated for hard- 
sphere glasses in Ref. $2(^. 
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The two-step crystallization process in hard-sphere glasses 
has been studied by Monte Carlo simulations by Schilling et 
al. II 1811 . at high volume fraction $ = 0.54. They found 
a rather gradual development of crystallinity from the disor- 
dered state as the simulation goes on, while in our case crys- 
tallization is rather sudden. 




FIG. 4: al) Snapshot from a simulation at ko = 3 in which the 
aggregate is still in the disordered state. Two colloids surrounded 
by a NaCl-type short-range order are shown as big spheres, together 
with their six nearest neighbors. a2) The final structure in the same 
simulation, showing a well-ordered NaCl-type crystallite, bl) Snap- 
shot from a simulation at na — 5 in which the aggregate is still in 
the disordered state. One colloid surrounded by a CsCl-type short- 
range order is shown as a big sphere, together with its eight nearest 
neighbors, bl) The final structure in the same simulation, showing a 
well-ordered CsCl-type crystallite. 



In the disordered phase, for na < 3.3, the range of the re- 
pulsion between colloids with the same charge is sufficiently 
long to prevent the occurrence of configurations in which 
more than six colloids of the same sign can surround one col- 
loid of the opposite sign. These configurations would be in- 
deed energetically favorable, because the CsCl phase is still 
the most stable from this point of view, but they are very dif- 
ficult to reach in an aggregation process starting from a di- 
lute volume fraction. This is the cause of the formation of 
crystallites of the metastable NaCl phase. Once formed, these 
crystallites can survive for long lifetimes. 

Rather surprisingly, we note that the formation of NaCl- 
type crystallites is indeed easier when this bulk phase is 
metastable. We have in fact tried aggregation simulations for 
Ka = 1 and U* = 9, a situation in which the NaCl bulk phase 
is more stable than the CsCl phase due to entropic effects JH . 
However, the long range of the repulsion between particles 
with the same charge renders the aggregation itself quite dif- 
ficult, so that we were not able to observe any crystallization 



within t ~ 210 3 t*. 

Finally we discuss the importance of hydrodynamic effects 
by comparing the results of BD and SRD-MD. 

The SRD-MD is much heavier computationally than BD so 
that it is very cumbersome to simulate the entire process of 
two-step aggregation and reorganization, because the nucle- 
ation of the ordered phase from the liquid can occur after very 
long times. For this reason we have decided to restrict the 
comparison to the analysis of the crucial parts of the process. 
First of all, we have started both BD and SRD-MD from initial 
randomly dispersed non-overlapping configurations and ana- 
lyzed the initial aggregation stage, leading to the formation of 
the disordered liquid-like phase, to check whether this phase 
forms on the same time scale and presents the same properties 
in both BD and SRD-MD. Then we have analyzed the for- 
mation of the ordered crystal, starting from configurations in 
which a NaCl seed is already well developed. These config- 
urations have been taken from previous BD simulations. The 
configurations with NaCl seed have been let evolve by BD and 
by SRD-MD. The results for na = 3, averaged over 4 simula- 
tions in each case, are reported in Fig. [5] where the behavior 
of Pi versus time is shown. The figure shows that both BD 
and SRD-MD produce the same kind of behavior. The liquid- 
like disordered phase is formed on comparable time scales in 
both cases, and it presents the same structure, with P2 oscillat- 
ing around 0.035-0.04. Also the growth of the NaCl seeds in 
the liquid occurs on the same time scale (in SRD-MD it seems 
only somewhat slower than in BD), producing the same kind 
of final structures. These results indicate that hydrodynamic 
effects do not play a major role in the two-step formation pro- 
cess of the ordered NaCl crystallites. 



CONCLUSIONS 

In this work we have analyzed the heteroaggregation of op- 
positely charged colloids by means of a combination of dif- 
ferent simulation techniques, including Brownian Dynamics, 
global optimization searches and Stochastic Rotation Dynam- 
ics - Molecular Dynamics. Aggregation in dilute systems has 
been considered. Our results show that colloidal aggrega- 
tion can be dominated by kinetic effects which can lead to 
the formation of very well ordered metastable crystals. These 
crystals significantly differ from the disordered metastable fee 
phases that were previously found in Ref. 11211 at higher vol- 
ume fractions. 

Our simulations show that, if the range of the potential is 
sufficiently long (na < 3.3), kinetic effects cause the forma- 
tion of crystallites of the metastable NaCl-type bulk phase, 
which thus can be obtained also for binary colloids of same 
size and balanced charge. These crystallites form within a 
two-step mechanism, in which the first step is the aggregation 
of a dense metastable liquid whose short-range order has pre- 
vailing NaCl features. This triggers the transition to ordered 
NaCl-type crystallites. Therefore the most likely outcome of 
the aggregation process is not the crystal structure of greater 
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FIG. 5: Comparison between BD (circles) and SRD-MD (triangles) 
results. Upper panel: behavior of P2 in the initial stages of the ag- 
gregation process, in which a liquidlike aggregate is formed starting 
from a diluted randomly dispersed configuration. Lower panel: be- 
havior of P2 for the growth of the NaCl crystal from a liquid-like 
aggregate in which a well-developed NaCl seed is already present. 
Time is measured in units of t* . 

thermodynamic stability, but the structure whose short-range 
order is closer to the short-range order in the metastable liq- 
uid. This indicates a mechanism of general character which 
is likely to be relevant for a variety of nucleation processes. 
The analysis of the most stable structures of small clusters 
has also shown that the formation of the NaCl phase is not 
related to the structures of very small aggregates, that are ei- 
ther amorphous or of CsCl structure. Our results show also 
that NaCl-type crystallites more easily form when the NaCl 
bulk phase is metastable or only marginally stable. Finally, 
the inclusion of hydrodynamic effects in the simulations has 
not produced major effects on the aggregation process. All 
these findings are important for designing strategies to obtain 
the desired structures of colloidal crystals in practical cases. 
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dio dell'eteroaggregazione di particelle colloidali ceramiche 
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